!! 基于蛙跳法求解时间步迭代，蛙跳法源于速度 verlet 法，其精度略低于 4 阶。

!> author: 左志华
!> date: 2022-06-10
!> version: alpha
!>
!> Leap-frog method for solving the Navier-Stokes equation <br>
!> 对实粒子进行蛙跳积分法
module sph_leap_frog

    use sph_kinds, only: rk
    use sph_region, only: region_t
    use sph_navier_stokes_solver, only: NS_fcn
    use sph_density_summation, only: density_summation_type
    use sph_env, only: vmin, rmin, umin, nums

    abstract interface
        !> Time integration <br>
        !> 时间积分方案
        subroutine time_integration_fcn(region, istep, t, dt, func)
            import :: region_t, rk
            type(region_t), intent(inout) :: region     !! region of particles <br>
                                                        !! 粒子宏观变量域
            integer, intent(in) :: istep                !! step number <br>
                                                        !! 步数
            real(rk), intent(in) :: t                   !! current time <br>
                                                        !! 当前时间
            real(rk), intent(in) :: dt                  !! time step <br>
                                                        !! 时间步长
            procedure(NS_fcn) :: func
        end subroutine time_integration_fcn
    end interface

contains

    !> Calculate the position, velocity, and acceleration of the first step <br>
    !> 计算第 2 步的位置、速度、加速度
    !> @note 本例程依赖预测半步长所需要的临时变量，如 rmin，umin，vmin。
    subroutine leap_frog_init(region, istep, t, dt, func)
        type(region_t), intent(inout) :: region     !! particle macroscopic variables domain <br>
                                                    !! 粒子宏观变量域
        integer, intent(in) :: istep                !! current step <br>
                                                    !! 当前步数
        real(rk), intent(in) :: t                   !! Current time <br>
                                                    !! 当前时间
        real(rk), intent(in) :: dt                  !! time step <br>
                                                    !! 时间步长
        procedure(NS_fcn) :: func

        call func(region, istep, t)                 ! 计算第 1 步/ 0 时刻的力 (加速度)，@todo 直接保存加速度量
        associate (n => nums(1), halfdt => dt/2)

            ! 精确量，物理量不同步
            region%particles%vel(:, :n) = region%particles%vel(:, :n) + region%particles%acc(:, :n)*halfdt      ! 速度前进半步长
            region%particles%u(:n) = region%particles%u(:n) + region%particles%du(:n)*halfdt                    ! 内能前进半步长
            if (density_summation_type == 1) then
                region%particles%rho(:n) = region%particles%rho(:n) + region%particles%drho(:n)*halfdt          ! 密度前进半步长
            end if
            region%particles%loc(:, :n) = region%particles%loc(:, :n) + region%particles%vel(:, :n)*dt        ! 位置前进一步长

            ! @todo 保存上时刻量，以便程序中断后可继续
            ! 预测量，使得物理量同步
            ! 保存半步长量
            vmin(:, :n) = region%particles%vel(:, :n)
            umin(:n) = region%particles%u(:n)
            if (density_summation_type == 1) then
                rmin(:n) = region%particles%rho(:n)
            end if

            ! 预测半步长，使得物理量计算力时，步长一致
            region%particles%vel(:, :n) = region%particles%vel(:, :n) + region%particles%acc(:, :n)*halfdt
            region%particles%u(:n) = region%particles%u(:n) + region%particles%du(:n)*halfdt
            if (density_summation_type == 1) then
                region%particles%rho(:n) = region%particles%rho(:n) + region%particles%drho(:n)*halfdt
            end if

        end associate
        call func(region, istep, t + dt)    ! 计算第 2 步的加速度

    end subroutine leap_frog_init

    !> Calculate the position, velocity, and acceleration of the second step and later <br>
    !> 计算第 3 步及以后的位置、速度、加速度
    !> @note 本例程依赖预测半步长所需要的临时变量，如 rmin，umin，vmin。
    subroutine leap_frog_calc(region, istep, t, dt, func)
        type(region_t), intent(inout) :: region     !! particle macroscopic variables domain <br>
                                                    !! 粒子宏观变量域
        integer, intent(in) :: istep                !! current step <br>
                                                    !! 当前步数
        real(rk), intent(in) :: t                   !! Current time <br>
                                                    !! 当前时间
        real(rk), intent(in) :: dt                  !! time step <br>
                                                    !! 时间步长
        procedure(NS_fcn) :: func

        associate (n => nums(1), halfdt => dt/2)

            ! 精确量，物理量不同步
            region%particles%vel(:, :n) = vmin(:, :n) + region%particles%acc(:, :n)*dt
            region%particles%u(:n) = umin(:n) + region%particles%du(:n)*dt
            if (density_summation_type == 1) then
                region%particles%rho(:n) = rmin(:n) + region%particles%drho(:n)*dt
            end if
            region%particles%loc(:, :n) = region%particles%loc(:, :n) + region%particles%vel(:, :n)*dt

            ! 预测量，使得物理量同步
            ! 保存半步长量
            vmin(:, :n) = region%particles%vel(:, :n)
            umin(:n) = region%particles%u(:n)
            if (density_summation_type == 1) then
                rmin(:n) = region%particles%rho(:n)
            end if

            ! 预测半步长，使得物理量计算力时，步长一致
            region%particles%vel(:, :n) = region%particles%vel(:, :n) + region%particles%acc(:, :n)*halfdt
            region%particles%u(:n) = region%particles%u(:n) + region%particles%du(:n)*halfdt
            if (density_summation_type == 1) then
                region%particles%rho(:n) = region%particles%rho(:n) + region%particles%drho(:n)*halfdt
            end if

        end associate
        call func(region, istep, t)             ! 计算第 n 步的力 (加速度)

    end subroutine leap_frog_calc

end module sph_leap_frog
